 use "${tmp}/Regfile_nafZE", replace

	
	xtset gp year
	
	
	destring ZE2010, replace
	
	// define sample 
	egen mlabexid = group(mainlabex)
	qui su mlabexid
	replace mlabexid = r(max) + ZE2010 if mi(mlabexid)
	

	reghdfe ln_eng_sbrut E_post11  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(ZE2010#year  gp)  cl(mlabexid)	
	bys gp: egen insampeng = max(e(sample))
	keep if insampeng
	 
	
	// Figure B3a	
	reghdfe ln_eng_sbrut HEtop_???? HEmid_????  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(gp ZE2010#year)  cl(mlabexid)	
	
	
	regsave using "${tmp}/tmp", replace
	
preserve
	use "${tmp}/tmp", clear
		keep if substr(var, 1, 1) == "H"
		split var, p(_)
		destring var2, replace force
		rename var2 year
		local N = _N+2
		set obs `N'
		replace year = ${NORMYEAR} if year == .
		replace coef = 0 if coef == .
		replace stderr = 0 if stderr == .
		replace N = N[1] if N == .
		replace var1 = "HEmid" in `=_N'
		replace var1 = "HEtop" in `=_N-1'
		
	
		gen ub95 = coef + invttail(N, 0.025)*stder
		gen lb95 = coef - invttail(N, 0.025)*stder
		sort year
		replace year = year + 0.1 if var1 == "HEtop"
		twoway (rcap ub95 lb95 year if var1 == "HEtop"  , lpat(solid) yline(0) xline(2010.5)) (scatter coef year, col(navy))  (rcap ub95 lb95 year if var1 == "HEmid"  , lpat(solid) yline(0) xline(2010.5)), legend(label(1 "Top tercile") label(3 "Middle tercile") ring(0) pos(11) order(1 3)) yscale(fextend) xscale(fextend) ///
		xlab(2004(1)2019, angle(40) nogrid) ylab(#6, nogrid) ytitle("Coefficient") xtitle("") ysize(9) xsize(16)
		graph export  "${outpath}/FigureB3a.pdf", as(pdf) replace
		
restore

	// Figure B3b	

	reghdfe ln_eng_sbrut HEtop_???? HEmid_????  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(gp ZE2010#year ape_n#year)  cl(mlabexid)	
	
	regsave using "${tmp}/tmp", replace
	
preserve
	use "${tmp}/tmp", clear
		keep if substr(var, 1, 1) == "H"
		split var, p(_)
		destring var2, replace force
		rename var2 year
		local N = _N+2
		set obs `N'
		replace year = ${NORMYEAR} if year == .
		replace coef = 0 if coef == .
		replace stderr = 0 if stderr == .
		replace N = N[1] if N == .
		replace var1 = "HEmid" in `=_N'
		replace var1 = "HEtop" in `=_N-1'
		
	
		gen ub95 = coef + invttail(N, 0.025)*stder
		gen lb95 = coef - invttail(N, 0.025)*stder
		sort year
		replace year = year + 0.1 if var1 == "HEtop"
		twoway (rcap ub95 lb95 year if var1 == "HEtop"  , lpat(solid) yline(0) xline(2010.5)) (scatter coef year, col(navy))  (rcap ub95 lb95 year if var1 == "HEmid"  , lpat(solid) yline(0) xline(2010.5)), legend(label(1 "Top tercile") label(3 "Middle tercile") ring(0) pos(11) order(1 3)) yscale(fextend) xscale(fextend) ///
		xlab(2004(1)2019, angle(40) nogrid) ylab(#6, nogrid) ytitle("Coefficient") xtitle("") ysize(9) xsize(16)
		graph export  "${outpath}/FigureB3b.pdf", as(pdf) replace
		
restore

	
	// Figure B3c	
	reghdfe ln_eng_hrs HEtop_???? HEmid_????  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(gp ZE2010#year)  cl(mlabexid)	
	
	
	regsave using "${tmp}/tmp", replace
	
preserve
	use "${tmp}/tmp", clear
		keep if substr(var, 1, 1) == "H"
		split var, p(_)
		destring var2, replace force
		rename var2 year
		local N = _N+2
		set obs `N'
		replace year = ${NORMYEAR} if year == .
		replace coef = 0 if coef == .
		replace stderr = 0 if stderr == .
		replace N = N[1] if N == .
		replace var1 = "HEmid" in `=_N'
		replace var1 = "HEtop" in `=_N-1'
		
	
		gen ub95 = coef + invttail(N, 0.025)*stder
		gen lb95 = coef - invttail(N, 0.025)*stder
		sort year
		replace year = year + 0.1 if var1 == "HEtop"
		twoway (rcap ub95 lb95 year if var1 == "HEtop"  , lpat(solid) yline(0) xline(2010.5)) (scatter coef year, col(navy))  (rcap ub95 lb95 year if var1 == "HEmid"  , lpat(solid) yline(0) xline(2010.5)), legend(label(1 "Top tercile") label(3 "Middle tercile") ring(0) pos(11) order(1 3)) yscale(fextend) xscale(fextend) ///
		xlab(2004(1)2019, angle(40) nogrid) ylab(#6, nogrid) ytitle("Coefficient") xtitle("") ysize(9) xsize(16)
		graph export  "${outpath}/FigureB3c.pdf", as(pdf) replace
		
restore

	// Figure B3d	

	reghdfe ln_eng_hrs HEtop_???? HEmid_????  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(gp ZE2010#year ape_n#year)  cl(mlabexid)	
	
	regsave using "${tmp}/tmp", replace
	
preserve
	use "${tmp}/tmp", clear
		keep if substr(var, 1, 1) == "H"
		split var, p(_)
		destring var2, replace force
		rename var2 year
		local N = _N+2
		set obs `N'
		replace year = ${NORMYEAR} if year == .
		replace coef = 0 if coef == .
		replace stderr = 0 if stderr == .
		replace N = N[1] if N == .
		replace var1 = "HEmid" in `=_N'
		replace var1 = "HEtop" in `=_N-1'
		
	
		gen ub95 = coef + invttail(N, 0.025)*stder
		gen lb95 = coef - invttail(N, 0.025)*stder
		sort year
		replace year = year + 0.1 if var1 == "HEtop"
		twoway (rcap ub95 lb95 year if var1 == "HEtop"  , lpat(solid) yline(0) xline(2010.5)) (scatter coef year, col(navy))  (rcap ub95 lb95 year if var1 == "HEmid"  , lpat(solid) yline(0) xline(2010.5)), legend(label(1 "Top tercile") label(3 "Middle tercile") ring(0) pos(11) order(1 3)) yscale(fextend) xscale(fextend) ///
		xlab(2004(1)2019, angle(40) nogrid) ylab(#6, nogrid) ytitle("Coefficient") xtitle("") ysize(9) xsize(16)
		graph export  "${outpath}/FigureB3d.pdf", as(pdf) replace
		
restore
	
	// Figure B3e
	reghdfe ln_eng_hrwg HEtop_???? HEmid_????  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(gp ZE2010#year)  cl(mlabexid)	
	
	
	regsave using "${tmp}/tmp", replace
	
preserve
	use "${tmp}/tmp", clear
		keep if substr(var, 1, 1) == "H"
		split var, p(_)
		destring var2, replace force
		rename var2 year
		local N = _N+2
		set obs `N'
		replace year = ${NORMYEAR} if year == .
		replace coef = 0 if coef == .
		replace stderr = 0 if stderr == .
		replace N = N[1] if N == .
		replace var1 = "HEmid" in `=_N'
		replace var1 = "HEtop" in `=_N-1'
		
	
		gen ub95 = coef + invttail(N, 0.025)*stder
		gen lb95 = coef - invttail(N, 0.025)*stder
		sort year
		replace year = year + 0.1 if var1 == "HEtop"
		twoway (rcap ub95 lb95 year if var1 == "HEtop"  , lpat(solid) yline(0) xline(2010.5)) (scatter coef year, col(navy))  (rcap ub95 lb95 year if var1 == "HEmid"  , lpat(solid) yline(0) xline(2010.5)), legend(label(1 "Top tercile") label(3 "Middle tercile") ring(0) pos(11) order(1 3)) yscale(fextend) xscale(fextend) ///
		xlab(2004(1)2019, angle(40) nogrid) ylab(#6, nogrid) ytitle("Coefficient") xtitle("") ysize(9) xsize(16)
		graph export  "${outpath}/FigureB3e.pdf", as(pdf) replace
		
restore

	// Figure B3f

	reghdfe ln_eng_hrwg HEtop_???? HEmid_????  if NAF_labex_accept == 1 & ZE_labex_accept == 1 , a(gp ZE2010#year ape_n#year)  cl(mlabexid)	
	
	regsave using "${tmp}/tmp", replace
	
preserve
	use "${tmp}/tmp", clear
		keep if substr(var, 1, 1) == "H"
		split var, p(_)
		destring var2, replace force
		rename var2 year
		local N = _N+2
		set obs `N'
		replace year = ${NORMYEAR} if year == .
		replace coef = 0 if coef == .
		replace stderr = 0 if stderr == .
		replace N = N[1] if N == .
		replace var1 = "HEmid" in `=_N'
		replace var1 = "HEtop" in `=_N-1'
		
	
		gen ub95 = coef + invttail(N, 0.025)*stder
		gen lb95 = coef - invttail(N, 0.025)*stder
		sort year
		replace year = year + 0.1 if var1 == "HEtop"
		twoway (rcap ub95 lb95 year if var1 == "HEtop"  , lpat(solid) yline(0) xline(2010.5)) (scatter coef year, col(navy))  (rcap ub95 lb95 year if var1 == "HEmid"  , lpat(solid) yline(0) xline(2010.5)), legend(label(1 "Top tercile") label(3 "Middle tercile") ring(0) pos(11) order(1 3)) yscale(fextend) xscale(fextend) ///
		xlab(2004(1)2019, angle(40) nogrid) ylab(#6, nogrid) ytitle("Coefficient") xtitle("") ysize(9) xsize(16)
		graph export  "${outpath}/FigureB3f.pdf", as(pdf) replace
		
restore